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We study imbalanced fermionic superfluids in an array of one dimensional tubes at the incipient 
dimensional crossover regime, wherein particles can tunnel between neighboring tubes. In addition 
to single-particle tunneling (ST), we consider pair tunneling (PT) that incorporates the interaction 
effect during the tunneling process. We find that with an increase of PT strength, a system of low 
global polarization evolves from a structure with a central Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) 
state to one with a central BCS-like fully-paired state. For the case of high global polarization, the 
center region exhibits pairing zeros embedded in a fully-paired order. In both cases, PT enhances 
the pairing gap, suppresses the FFLO order and leads to spatial separation of fully-paired and fully- 
polarized regions, same as in higher dimensions. Thus, we show that PT rather than ST is essential 
for the development of signatures characteristic of the incipience of the dimensional crossover. 

PACS numbers: 37.10.Jk, 67.85. Lm, 71.10.Pm 



Superconductivity and ferromagnetism are two ubiq- 
uitous but competing phenomena in condensed matter 
systems. Spin imbalance and magnetic field induced by 
ferromagnetism tend to suppress Cooper pairing, which 
is responsible for superconductivity. For more than four 
decades an interesting phase, the Fulde-Ferrell-Larkin- 
Ovchinnikov (FFLO) state [H [2] , has been suggested as 
the concurrence of both ferromagnetic and oscillatory 
superconducting orders [21 H], but its direct confirma- 
tion is still elusive. Recently, owing to the capability 
of controlling particle densities, tuning interactions and 
cooling into quantum degeneracy [5] E] , cold atomic sys- 
tems have become a promising platform for searching 
for FFLO order [5] . In experiments of trapped Fermi 
gases, density profiles that reflect the interplay of spin im- 
balance (ferromagnetic order) and Cooper pairing have 
been observed [5HTT] . In addition, experiments have also 
revealed significant dimensional dependence of the pro- 
files: in three dimensions (3D) a fully paired profile takes 
place at the trap center and a polarized profile does off 
center [HI [TO], while in one dimension (ID) the central 
region is always polarized [TT]. These observations agree 
with theoretical studies of the FFLO state in ID and 
the trap-induced phase separation in 3D [I2HI9] , but the 
marked difference between these two limits also raises the 
need for understanding the intermediate regime. Several 
works have focused on the dimensional crossover regime 
of various kinds of continuous systems [2"0H2"5] or Hub- 
bard lattices 26, 27J, but with a different emphasis than 
the present work. 

In this paper, we study a realizable system of a two 
dimensional optical lattice array of ID tubes, subject 
to a global trapping potential [TTJ [251 12S]- The incipi- 
ent dimensional crossover regime of this system, which 
can be experimentally accessed by lowering the lattice 
depth, is modeled by incorporating the kinetics of single- 
particle tunneling (ST) as well as a key ingredient rep- 




FIG. 1: (color online) Illustration of single particle and pair 
tunneling processes, ti and t2, respectively, between two 
neighboring tubes in the array. On each tube, two opposite- 
spin atoms can form a bound pair (circled in the graph) in 
the presence of an attractive interaction g. The circled pair in 
the middle indicates that the two atoms remain paired during 
the t2 process. 



resenting the tunneling of paired opposite-spin atoms - 
pair tunneling (PT) - between neighboring tubes (as il- 
lustrated in Fig. [IJ . The ST leads to an interesting mag- 
netic compressible-incompressible phase transition anal- 
ogous to that in the Bose-Hubbard model (discussed in 
Ref. [2S]) but is not responsible for certain observed 
signatures in the profiles at the dimensional crossover 
regime (that will be shown below). By considering PT, 
we are able to describe the evolution of profiles from ID 
toward 3D and obtain the emerging signatures of the di- 
mensional crossover at various global polarizations, such 
as the inversion of the fully-paired and polarized centers 
as well as the growing spatial separation between the 
fully-paired and fully-polarized regions. Following, we 
first discuss the microscopic physical cause of the PT and 
evaluate its strength using a two-channel model. We then 
construct a model Hamiltonian and apply a Bogoliubov- 
de Genne (BdG) treatment to solve for the density and 
pairing profiles of the system. Finally, we present and 
discuss the results. 

For free fermions in a lattice potential, the inter-site 
kinetics is well described by ST processes of strength 
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t\ [HUj. For an attractively interacting case where two 
opposite-spin atoms form a pair of binding energy e&, the 
kinetics of the pair tunneling would, in principle, be in- 
corporated as a process in which the two atoms of a pair 
split, separately tunnel to the other site and re-bind. This 
is contained in the ST model as a second order process of 
strength tf/eb and accounts for the Josephson phenom- 
ena in the presence of superffuid orders, cf. [22] . 

However, in cold-atom experiments, the interaction is 
induced via a Feshbach resonance [51] . which is con- 
trolled by the tuning of a magnetic field affecting the 
hyperfine energy splittings. Because the field is ap- 
plied throughout the system, the interaction that leads 
to pairing exists on and between lattice sites. There- 
fore, we expect that atoms that remain paired during 
the whole tunneling event can be another viable process 
(see illustration in Fig. [I]). Such process is described 
as tunneling of the paired atoms, with strength denoted 
as t 2 . One can estimate t 2 around the Feshbach res- 
onant regime using a two-channel model that incorpo- 
rates atomic and molecular degrees of freedom (see Ap- 
pendix) . By integrating out the molecular degrees of free- 
dom, we obtain ti — a 2 t m , where t m is the molecular 
tunneling and a is the ratio of the molecule-atom cou- 
pling strength to the molecular chemical potential. Con- 
sidering the tunneling strength in optical lattices given 
by -^E R (V /E R ) 3/4 eM-^VWE^] with V being the 
optical-lattice depth and E R being the recoil energy [32] , 
we find t 2 /ti — y/2a 2 exp[— 2^V /E R ]. This expres- 
sion shows that the PT effect becomes more important 
and can even overcome the ST as the lattice depth de- 
creases. In this same limit, perturbative-RG analysis 
starting from ID becomes unreliable. In addition, our 
system shows prevasive pairing effects and is thus differ- 
ent from an array of Josephson junctions that lack the 
pairing mechanism in the insulating barriers between the 
superconductors. Therefore, by providing extra channels 
for Josephson tunneling, the PT processes can lead to an 
enhancement of the superfluid order and its cross-tube 
coherence, anticipating the dimensional crossover regime. 

Incorporating both the ST and PT effects, the tube 
lattices occupied by up-spin (majority) and down-spin 
(minority) atoms, Vv=tA>r( 2 )i are hence described by 
the microscopic Hamiltonian 

Jz r a 

+ I E HlE^'^-^tr^lr'W t r) 5 (l) 
Jz ( rr <) a 

with the z direction along the tube's axis and r = (x, y) 
denoting tube indexes in the plane perpendicular to z. 
The one-particle Hamiltonian, H® — — {h 2 /2m)d 2 + 
m(ui 2 r 2 + ujz 2 )/2 ~ /j, a , includes the kinetic energy in 
the z direction, the global trapping potential and the 



spin-dependent chemical potentials. The on-tube atomic 
coupling, g (taken positive for attractive interaction), is 
for the ID-tube limit [33]. The ST (PT) of strength t 1 
(t 2 ) j takes place between nearest-neighbor tubes, (rr'). 

Applying the BdG mean-field theory [33], (which has 
successfully described tube lattices without PT [25] and 
a variety of tube confinements [HI H5 ] fT9 ] I2T ] 135 ] I3B ] ). 
we construct a mean-field Hamiltonian, Hm, by corre- 
spondingly replacing the quartic operators in Eq. |l]) with 
quadratic ones coupled to three different mean fields, 

H ™ = / E [ E tin: {K + U„) $„ 

<rr'),cr 

Here the Hartree field, U ar (z), and the BCS gap field, 
A. r (z), are standard variational fields in previous BdG 
studies. We introduce a tunneling field, T ar r' ( z ), as a new 
ingredient to describe the effective tunneling under the 
influence of both t\ and t 2 . We rotate Hm into the quasi- 
particle basis, 7„, through a Bogoliubov transformation, 
4><rr(z) = Y, n [ u nav{z)^ niJ - av* nrjr {z)^ nS \ (where a = 
—cr), and derive extended BdG equations for the quasi- 
particle wave functions, u n<7 and v na , as well as the corre- 
sponding energies, e na . The condition, S (H — TS) = 0, 
which guarantees solutions of an equilibrium state at 
temperature T, leads to the self-consistent relations, 



U„ = 



A r = 



Ttt' u 



-dypir^r) = -gYWriarffna, (3) 

n 

r' 

E ( -9UntrV* lT ~ *2 E W «tr'«Inr' ) JW>( 4 ) 



(5) 



-ti-t 2 {i>l r ,i 



where f n<7 — [exp(e n(J / 'k R T) + I] -1 is the Fermi distri- 
bution function and , runs over all tubes at r' next 
to r. Equation Q shows that the magnitude of the 
pairing gap is enhanced by t 2 in uniform lattices where 
(VvVv) = vPir'^fr') (and would also be in trapped sys- 
tems, as expected through a local density approximation 
argument). This enhancement tends to stabilize the fully 
paired phase against being invaded by unpaired majority 
atoms; in analogy to the Meissner effect [37] which pre- 
vents the superconducting bulk from being penetrated 
by the magnetic field. When t 2 — 0, T = —t\ turns 
Eq. ([2| back to that for the Hamiltonian with only ST 
(discussed in Ref. [S3]). We numerically solve the BdG 
equation and apply the solutions to calculate the spatial 
profiles of pairing gap A, total density p = p^ + p^, and 
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FIG. 2: (color online). Axial profiles (AP) of total density 
p, spin imbalance s, (solid light blue and dark red curves, 
respectively, axis on the left of graph), and average magni- 
tude of pairing gap |Aj (dashed green curve, axis on the right 
of graph) for various pair-tunneling strengths t2, and global 
polarizations P. Rows from top to bottom correspond to 
t2 = 0, 0.014ef,, and 0.05et, respectively; while the left and 
right columns correspond to P — 0.25 and 0.5, respectively. 
On the left column the dash-dotted lines demarcate regions of 
(i) FFLO, (ii) BCS-like fully-paired, and (iii) fully-polarized 
states. The data obtained were for systems of 2400 parti- 
cles in a 10 x 10 tube array with global trapping frequency 
uj — 0.0625ef,/ft (which defines the oscillator length, a), single 
particle tunneling ti = 0.014eb, and temperature T = O.lei,. 
These parameters are similar to those used in experiments. 
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FIG. 3: (color online). Value of the gap function (in units of 
hw) at the center of each tube (z — 0) for various t% and P 
(convention as presented in Fig. [21. Here we show data for 
5x5 tubes in the fourth quadrant of the 10 x 10 tube array, 
in which the top left entry of each panel corresponds to the 
most central tube. The other quadrants are similar due to 
4-fold rotational symmetry. 



spin imbalance (or magnetization) s = p^ — where 
p a = (V4.-0O-) is the density profile of a species. Below 
we present the results for a spherically trapped system 
(uj r = us). 

We look at the influence of t 2 at fixed t\ — 0.014eb, the 
latter corresponding to a typical lattice depth of 7 En and 
thus into the dimensional crossover regime. In Fig. [2] we 
plot the axial profiles of p, s and the average of |A| by 
tracing out the r degree of freedom. The first and sec- 
ond columns correspond to a lower global polarization of 
P = 25% (LP) and a higher one of P = 50% (HP), re- 
spectively [Mj. From top to bottom rows, is chosen to 
be either zero, comparable to t%, or large compared with 
t\, respectively. We see that in the LP case at t% = 0, 
the axial profile exhibits (i) an FFLO center with oscil- 
latory A, (ii) a BCS-like shoulder with non-oscillatory 
A, and (iii) a normal tail having zero A. At the in- 
termediate t 2 value, this tri-layered structure remains. 
However, the FFLO center shrinks, the BCS-like region 
extends toward the center accompanied with a drop in 
imbalance, and the normal tail grows. This indicates a 
transfer of unpaired majority atoms from the center to 



the tail. We notice that the gap profile develops small 
ripples between the BCS-like shoulder [(ii)] and the nor- 
mal tail [(iii)], suggesting the incipience of an FFLO layer 
[(i)] here. At the large t 2 value, the FFLO center is com- 
pletely conquered by the BCS-like state and disappears, 
leaving a large fully-polarized tail and a thin FFLO layer 
in between them. Because the FFLO and BCS centers 
are distinctive of ID [11] and 3D [9j[T0] trapped systems, 
respectively, this result shows the evolution of the sys- 
tem from ID toward 3D, driven by t 2 (compared with 
increasing ti). 

In the HP case, the system always has a center with 
oscillatory A and a fully-polarized normal tail. In the 
oscillatory-pairing region, the imbalance profile exhibits 
characteristic out-of-phase oscillations, with the concur- 
rence of local minima (maxima) of s and local maxima 
(minima) of |A|. This behavior is due to the competi- 
tion between supcrfluid and ferromagnetic orders. An 
increase in t 2 enhances this competition, augmenting the 
magnitude of the out-of-phase oscillations and repelling 
a portion of the unpaired majority to the normal tail re- 
gion. At large t 2 (— 0.056;,), the oscillations are large 
enough that the minima of s are almost zero. Such case 
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FIG. 4: (color online) The combined fraction of particles in 
the highly paired (s/p < 5%) and highly polarized (s/p > 
95%) regions, 7, vs t-z (left panel) or t\ (right panel) with 
the other fixed. The solid blue, dashed red, and dotted green 
curves represent cases with global polarization P = f2.5%, 
25% and 50%, respectively. 



is less like an FFLO state (oscillatory pairing accompa- 
nied with finite polarization), but more like spatial al- 
ternation of fully-paired superfluid and highly polarized 
normal gas. This phenomena, analogous to the phase 
separation in the LP case, is taken as a signature of the 
dimensional crossover between ID and 3D at higher po- 
larizations. We also notice that the structure of the pro- 
files is reminiscent of that of a system with vortex cores 
embedded in a superfluid bulk. 

We find that PT affects the pairing order not only 
along but also across the tubes. Figure [3] shows the value 
of the gap function at the center of each tube (z = 0) 
in a 10 x 10 tube array. The left (right) column corre- 
sponds to the LP (HP) case, while, from top to bottom, 
rows correspond to zero, intermediate, and large t-i, re- 
spectively (as in Fig. |2|. We see that at zero £2 the sign 
of A changes, indicating an oscillatory behavior across 
the tubes. In the LP case when increases, the oscil- 
lating nodes appear in a more off-center region, as dis- 
cussed for the axial profiles along the tubes. At large t^, 
there is no oscillation of A across tubes in both LP and 
HP cases, showing the suppression of FFLO order. We 
notice in Figs. [2] and [3] that t-i enhances the maximum 
magnitude of the gap function, as expected from Eq. Q. 
This enhancement raises the critical temperature above 
which the pairing order vanishes and hence agrees with 
the increase of the superfluid transition temperature in 
quasi- ID systems |22j . 

Finally, we look at the phase separation of fully- 
paired and fully-polarized regions as a function of 
t%. We consider the combined fraction of particles 
in the highly paired {s/p < 5%) and highly polar- 
ized (s/p > 95%) regions of the axial profiles; 7 = 
f z p[6(0.05- s/p) + 0(s/p- 0.95)]/ f z p, where 9 is the 
step function. The larger 7 is, the stronger the phase 
separation the system shows. The left panel in Fig. [4] 
shows that 7 monotonically increases with ti at three 
various polarizations when t\ is fixed. In the cases of 
P = 12.5% and 25% the sudden changes indicate the oc- 
currence of the BCS-like center replacing the FFLO cen- 
ter. For comparison we plot also 7 vs t\ at fixed ti in the 



right panel and observe that 7 shows almost no change at 
the three polarizations. This result highlights that it is 
£2, rather than tx, that accounts for the phase separation 
and hence is essential for the correct modeling describing 
the physics at the incipience of the dimensional crossover 
regime. 

In conclusion, considering the microscopic physics of 
cold atomic systems, we have incorporated both ST and 
PT processes to effectively model imbalanced fermionic 
superfluids in an array of ID tubes at the incipience of 
the dimensional crossover. Our calculations show that 
the PT strength is a main factor for the evolution of the 
system profiles deviating from the ID limit, which ex- 
hibits a central FFLO state, toward the development of 
3D signatures, including a central fully-paired state in the 
LP case and spatial separation between fully-paired and 
fully-polarized states in both LP and HP cases. These 
features are reflected in the directly observed density pro- 
files and the pairing orders which can be probed in time- 
of-flight experiments ^38-41 . Our model can be easily 
generalized to incorporate higher-order, higher-hand or 
inter-band processes j3H 23] , which are expected to be of 
further help to investigate the system's transition to the 
continuous 3D limit. 
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APPENDIX: EVALUATION OF THE PAIR 
TUNNELING 

In this section we evaluate the PT strength using a 
two-channel Feshbach resonance model [6] [31] . Interact- 
ing fermionic atoms near a Feshbach resonance shall be 
described using both atomic and molecular degrees of 
freedom, tp a and <fi, respectively. In optical lattices, the 
partition function of the system is 

Z= [ D{^,^ ai }D{ ( j ) . l ,^}e-f dT ^ +s -\ (6) 



where S a contains terms associated only with the atomic 
degrees of freedom and 

S m = —tm (j>i<pj - fi m 4>i(f>i 

(ij),a i,(7 

+U am {falPiiM^i + H.c.) (7) 
i 

involves the molecular tunneling i m , molecular chemical 
potential /i m , and the atom-molecule coupling U am . Here 
we assume the molecular interaction is weak such that 
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the mean-field approximation, 4>i<j)i<f)i<t>i — > \4>i4>i) 4>i4n, 
can be applied to incorporate the interaction as effective 
contributions to the chemical potentials. We integrate 
out the molecular variable <$> in Eq. ([6]) and obtain 



(8) 



where S' a is expanded as 



s: = 



ir 



tn 
fir, 



- 2^ ^i^iii^iii'^ 

i 



(9) 



The first term in Eq. ^ can be treated as an effective 
inter-atomic interaction, while the second one appears as 
PT. Therefore, we obtain 



U, 



2 

am , 



(10) 



or a = U am / \i m as mentioned before in the text. From 



Eq. ( 10 1 we see that even if the molecular tunneling is 
smaller than the atomic tunneling (t m < t\) 1 ti can be 
comparable with or even larger than t\ in some param- 
eter regimes. Notice that we focus here on the leading 
processes in the lowest Bloch band. Incorporating higher- 
order, higher-hand or inter-band processes is one of the 
immediate extensions of our model 1421 EES- 
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